Signed Simplicial Decomposition and Overlay of 
n-D Polytope Complexes 

Norbert Paul 



Abstract Polytope complexes are the generalisation of polygon meshes in geo- 
information systems (GIS) to arbitrary dimension, and a natural concept for ac- 
cessing spatio-temporal information. Complexes of each dimension have a straight- 
forward dimension-independent database representation called Relational Complex. 
Accordingly, complex overlay is the corresponding generalisation of map overlay in 
GIS to arbitrary dimension. Such overlay can be computed by partitioning the cells 
into simplices, intersecting these and finally combine their intersections into the re- 
sulting overlay complex. Simplex partitioning, however, can expensive in dimension 
higher than 2. In the case of polytope complex overlay signed simplicial decompo- 
sition is an alternative. This paper presents a purely combinatoric polytope complex 
decomposition which ignores geometry. In particular, this method is also a decom- 
position method for non-convex polytopes. Geometric n-D-simplex intersection is 
then done by a simplified active-set-method — a well-known numerical optimisation 
method. "Summing" up the simplex intersections then yields the desired overlay 
complex. 
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1 Introduction 

An important query operation in GIS is the overlay of some given "topologies" to 
generate new such "topologies". An example could be cadastral land-owner data 
overlaid with environmental stress data which helps to identify which owner is af- 
fected by what averse environmental influences. This operation seems to be missing 
in 3D cadastral applications: 
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However, 3D data management and analysis such as querying, manipulation, 3D map over- 
lay, 3D buffering have been largely neglected in spatial database systems and Geographic 
Information Systems. 1171 

To compute such overlay by first triangulate the input, then overlay these triangles 
and finally recombine the resulting intersections into the desired overlay complex is 
possible in 2-d, but such triangulation in 3-d, however, has Q (n 2 ) space complexity 
in general [7| and, hence, is very expensive. 

The special case of convex 3-d-shapes, however, is almost trivial: Fix an arbitrary 
interior point and make it the apex of a family of cones atop the boundary faces. If 
these faces have been triangulated before the result is a simplicial decomposition of 
the convex shape. 

Also the area of an arbitrary non-convex planar polygon can be computed by 
a signed sum of the triangle areas made up of one edge and a fixed arbitrary ver- 
tex in the polygon plane. These observations will here be generalised to arbitrary 
dimension and used to compute n-d complex intersection. 



2 Related Work 

Much work has been done on polytope decomposition, volume computation, and 
intersection — mostly, however, on convex polytopes represented by vertices (via 
convex hull) or by (intersecting) half spaces. That representations do not allow non- 
convex polytopes and some of the above mentioned problems (even vertex enumera- 
tion) are NP-hard if no fixed dimension upper bound is given [ 1 3 1 . By using polytope 
complexes instead, vertices, edges, faces etc. are already explicitly enumerated as a 
precondition and the above problems can be avoided. 

Volume computation of convex polytopes is discussed in ifPfll and in |4|where the 
latter also introduces signed simplicial decomposition. 

Most work on decomposition is dedicated to unsigned decomposition as studied, 
for example, in the survey [9 1. Unsigned simplicial decomposition of convex poly- 
tope complexes are known as Boundary Triangulation [4], and as Cohen Hickey's 
Triangulation [ 10], cited by J3J. 

A work very similar to this article is described in [5 1 and in [6|as "alternate hi- 
erarchical decomposition" (AHD). This is also a singed decomposition into convex 
parts and its authors, too, consider it useful to compute intersection, union and sym- 
metric difference between non-convex polytopes. However, with the shape 
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that method may not terminate. 



3 Basic Notions and Data Model 

In many cases spatial data can be considered a model of some partitioning of the 
two- or three-dimensional space in spatial "chunks" like areas, volumes, faces etc. 
just like Computer-Aided Design (CAD), GIS, 3-d city models, or subsoil geology 
models do. 

If such partitioning undergoes changes in time this can be considered a parti- 
tioning of the four-dimensional space-time into what might then be called "hyper- 
chunks". A volume, for example, may extend over a time interval and then be split 
in two at a time point after which two such volumes start to exist. Then that volume 
at the interval before the split can be considered a four-dimensional "hypervolume" 
bounded by two volumes which mark the splitting event as shown in Figure Q] 




Fig. 1 The spatio-temporal process of a volume object (left) splitting into two (middle) and one of 
them slipping away over time (middle to right) as a 4-d space-time-complex. 
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3.1 Topological Extension of the Relational Model 



What follows is a brief introduction into our topological data model [3 1. Consider 
our aforementioned subdivision of space (or space-time) satisfying the following 
preconditions: 

• A compact subset of the n-dimensional real space IR" is subdivided into a finite 
number of parts. 

• Each part is a connected /-dimensional manifold without boundary. 

• Each such part is fiat. That means it is within an /-dimensional affine subspace 
ofIR". 

• The boundary of each /-dimensional part is the union of parts of lower dimen- 
sions. 

We call this subdivision a finite polytope complex. Then the natural topology of the 
underlying n-dimensional real space IR" generates a so-called quotient topology. As 
the number of parts is finite, their topology is finite, too, and, hence, it is a so-called 
Alexandrov-topology [1 1 which has an important characteristic: It can be stored by 
a relation called "incidence graph" and so it fits into a relational database: 

Definition 1 (Topological Data Type). Let X be a set and R CX xX a relation on 
X. We call the pair (X,R) a topological data type. A subset A C X is said to be open 
in (X,R), iff all x £ X and a £ A satisfy xRa => x £ A. The relation R is also called 
the incidence relation of (X,R). 

So far we have only relabelled what is commonly known as simple directed graph. 
But note that every topology for a finite set can be stored in that simple manner. 

Example 1 (Combinatorial Square). We partition a unit square [0,1] x [0,1] C 
M 2 into nine elements: four vertices a — {(0,1)}, b = {(1,1)}, c — {(0,0)}, and 
d = {(1,0)}, four edges e = ]0,1[ x {1}, / = {0} x ]0,1[, g = {1} x ]0,1[, and 
h = ]0, 1 [ X {0}, and the face F = ]0, 1 [ x ]0, 1 [, which gives the topological space 
depicted at the left-hand side of the diagram below. At its right-hand side there is the 
corresponding topological datatype (X,R) with point setX = {a,b,c,d,e,f,g,h,F} 
and incidence relation 

R = {(x,y) | The right image has an arrow x — > y.} . 




a -s e 
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The justification for relabelling "graph" to "topological data type" is our adaption 
of continuous maps to topological data types: 

Definition 2 (Continuous Database Map). Let (X,R) and (Y,S) be topological 
data types and / : X — >■ Y a map. S* denotes the transitive and reflexive closure 
of S. We then call / a continuous database map iff (f(a),f(b)) € S* holds for all 
(a,b) £R. Then we write/: (X,R) ->■ (Y,S). 

The continuous database maps are exactly the continuous maps between the corre- 
sponding Alexandrov spaces. What we have in mind is a relational database table X 
of spatial entities together with a table R as an n:m relation type from X to itself. 



3.2 Algebraic Topology in Relational Databases 

We now extend our above data model to algebraic topology: A Relational Chain 
Complex is a topological data type (X,R) if the relation R also carries additional 
information about the orientation of each cell by attaching signs to the cell-cell- 
incidences. 

Here our partitioning of the space consists of a sequence (C„, . . . ,Cq) of sets C, 
of /-dimensional manifolds which together make up the entire space and which we 
call by abuse of language "/-cells". Note that the above preconditions guarantee, 
that at least our 0- "cells" and our l-"cells" are cells indeed. Each i+ 1-cell c is 
orientable and bounded by a set D of /-cells. By fixing an orientation for each cell 
in our partitioning we specify a sign for each cell d in D which depends on c: It is 
positive if c is touched by d's "front" side, negative at its "rear" side, or zero if d 
touches c by "both" sides. We refer to [12, p. 233] for details on what "front" and 
"rear" in arbitrary dimension mean. 

The above mentioned sign defines a function 







: "front" 


a c :D^{-l,0,+l};d^ < 


!• 


: "both" 






: "rear" 



Now we extend our relation R from the topological data type (X,R) to a function 

M : R ->• Z, (xi+i ,x t ) i-> a Xi+i (xi) (2) 
and store it into a database table M with schema 

M ( cell :X, boundar y :X, sigma:Z). 



Considering M a sparse matrix the matrix product of M with itself can be computed 

by 
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create view M_squared as 
select Ml. cell, M2. boundary, 

sum(Ml. sigma * M2. sigma) as sigma 
from MM1, MM2 
where Ml. boundary = M2. cell 
group by Ml. cell, M2. boundary. 



If M_squared only contains zero entries it denotes a complex boundary. The above 
SQL-statement defines the multiplication of two sparse matrices Ml and M2 by sim- 
ply replacing the f rom-clause by 



from Ml, M2. 



Now we can define our data model: 

Definition 3 (Relational Complex). A sequence (X„,...,Xq) of finite sets — the 
cells — together with a sequence (D„,...,Dj) of sparse X; x -matrices — the 
boundaries — is called a relational complex of dimension n if every matrix product 
Z), + i • Di of two consecutive matrices has only entries of value zero. 

This definition fixes a static dimension upper bound n. Dynamic dimension is also 
possible: Collect all cells into a table X, all matrices into a table D, and specify 
dimension merely by an integer attribute. To specify geometry we simply store the 
n coordinates of each vertex, say, in case n — 4, by x, y, z, and t. 

Example 2 (Combinatorial Square). The following relational complex for Example 
Q]has Edges e and h running from left to right whereas Edges / and g point down- 
wards. Face F gets a counter-clockwise orientation: 

The right-hand side shows the corresponding relational complex with point sets 
Xq = {a,b,c,d}, X\ = {e,f,g,h}, andX2 = {F} and boundaries 

Ri= I (x,y, ±1) | xeXi,y e X/_ i , arrow x — — » y in right image below, I . 





Note that the orientation of Edge h is "compatible" with the orientation of Face F 
so the entry in the boundary table is the tuple (F,h, +1), whereas g runs contrary to 
the face orientation as entry (F,h, — 1) indicates. 
Now consider all paths from face F to vertex a: 
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F s~a and F *-f ^a. 

When we multiply the coefficients of the first path we get ( — 1 ) x ( — 1 ) = + 1 and for 
the second path we have (+1) x (— 1) = — 1. The sum of both is zero — the ,F, a-entry 
of the matrix product. 

We also need a mapping between two relational complexes: 

Definition 4 (Relational Complex Morphism). Let (X„,...,X Q ) with bound- 
aries (D n ,...,D\) be a relational complex, and let (Y„,...,Yq) with boundaries 
(B n ,...,B\) be another relational complex. Then we call a sequence (F„,...,Fq) 
of (sparse) X{ x Y,-matrices a relational complex morphism, if each D, = FjBi 
except for zero entries. 

The above used difference A — B of sparse matrices pads missing entries in one of 
the matrices with zero. From now on we say "algebraically equal" instead of "a = b 
except for zero entries". 



4 Complex Overlay 

Now having presented our data model we introduce the problem: 

Two polytope complexes C and K partitioning a compact part \C\ and \K\ of the 
real vector space into flat manifolds have a common refinement of the intersection 
|C| n \K\: the polytope complex of |C| (~l \K\ of the non-empty pair-wise intersections 
of cells in C and K. We call this intersection overlay and denote it by C n K. 

Additionally there is a common refinement of |C| U \K\ which takes into account 
the exterior of the a complex if needed. This is union overlay, denoted by C U K. 

Example 3 (Two Squares and a Triangle). The union-overlay and intersection- 
overlay of a complex S which consists of two meeting squares ABDE and BCEF 
with another complex T — simply the triangle XYZ — may be the following com- 
plexes S U T (left) and S D T (right): 



X 



A 




Z 



D E F 

The problem is: How can the intersection complex be computed? 
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Triangulation is the partitioning of polygons into triangles which practically is in 
O(nlogn), theoretically even in 0(n) J8). 

The 3-d-analog is partitioning a polyhedron into tetrahedra and is often called 
"tetrahedralisation". The worst case lower bound on the number of resulting tetra- 
hedra is in Q. (n 2 ), where n is the number of cells Q. 

The corresponding dimension-independent notion is "simplicial decomposition". 
We know this is expensive but, happily, signed simplicial decomposition is often an 
alternative. 

We will now present combinatorial n-d signed simplicial decompositions which 
generalise the well-known methods Boundary Triangulation [4 1, and Cohen Hickey's 
Triangulation ([10|, cited by [4|) by only using algebraic and topological informa- 
tion provided as relational complex. 



5.1 Signed Boundary Decomposition 

A signed decomposition of a cell c into simplices ±A\ , . . . , ^A„ is a linear combi- 
nation of simplices A\, . . . ,A n such that 

c = aiAi-\ Ya„A„ (3) 

where a, ■ — +1 if the simplex is added and a, ■ = — 1 if it is subtracted. As the or- 
der in which we add or subtract simplices does not matter one might also call this 
"Commutative Constructive Solid Geometry" (CCSG). 

The algorithms presented here operate on the relational representation of a com- 
plex but we now take its "classical" view where each C, is the free Abelian group of 
cells in X/, and the boundary operators <9, are linear maps between them satisfying 
<?,_i (di(x)) = for all x E Q. Now assume such complex 

^cAQ-^-^Q (4) 

be given. 

To triangulate an n-cell c„ from C„ (hence c n G X n ) we simply add an arbitrary 
point a from the interior of c„ as new vertex a ("apex") to the vertices and replace 
our cell c„ by a "cone" a®d{c n ) over its boundary with apex a. The tensor operator 
a ®x is simply the concatenation of a and x, either by creating a pair (a,x) or by 
appending tuples: 

a (g> d (c n ) — a (g) (ai • dy -\ \-a m -d m ) 

= a ® 05i • d\-\ Y a®a m - d m (5) 

= ai • a ® d\ -\ V a m -a®d m 
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The boundary of one such such element a <E> di can be defined as 

d(a®di):=di — a®d(di). (6) 

This is a boundary operator, indeed. It slightly modifies the Eilenberg-Zilber- 
formula ifTTl for tensor product boundaries 

8{a®x) :=8{a)®x+{-\f ima a®8{x) (7) 

to 

d{a®x) :=d(a)(g>x-(-l) dima a<E)d(x) (8) 

and specifies d(a) = () which means that the boundary of the apex a is the empty 
tuple — the identity element of concatenation: (} ®x =x =x&) (}. Subtraction instead 
of addition is achieved by simply "shifting" the dimension of a from zero to 1 and 
considering (} having dimension 0. 

In contrast to the original, our variant of the Eilenberg-Zilber-formulais compat- 
ible with the simplicial boundary operator which we denote here by 8: 

d(a<g> (b,c,d,...)) 

= (b,c,d, . . .) - (a,c,d, ...) + {a,b,d, ...) — (a,b,c,.. .) H 

= 8((a,b,c,d,...}) (9) 

This even works for arbitrary n-simplices (ao, . . . ,a n ) (g> (b,c,d, ...). 

We remind that the simplicial boundary 8 of a simplex (vq, . . . , v n ) is defined as 

n 

<5((v ,...,v„)) :=J^(-iy(v ,...,Vi-i, v i+l ,...,v„) . (10) 

i=0 

Our "cone" over the boundary in fact has the same algebraic boundary as the 
original cell: 

d(a<gid(c„)) = d{c„) —a®d{d{c„)) = <9(c„) -a®0 

= d(c n ). (11) 

So a property of c„ computed via its boundary can also be computed by the proposed 
triangulation — even when a is outside c„ or a vertex of the boundary of c n . 

Example 4 (Triangulating a Square). The left hand side below shows a square S in 

counter-clockwise orientation (as indicated by the bent arrow ■>- and the right 

hand side its triangulation with a newly introduced apex "0": 
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The signs of the four triangles correspond to their orientation: If a triangle has clock- 
wise orientation (like (0,1,2) which visits its vertices 0, 1, and 2 in clock- wise 
order) it gets a negative sign. Each common edge of two adjacent signed triangles 
vanish when the boundary of their sum is computed. In <3((0, 1,3) — (0,1,2)) = 
+ (1,3) - (0,3) + (0,2) - (1,2) +0 x (0,1), for example, the common edge (0,1) 
cancels out. 

The boundary of the left hand side complex is 



d 2 (S) = -(l,2) + (l,3)-(2,4) + (3,4) 



(12) 



for the face and 



<M(1,2)) = (2)-(1) 
a«2,4» = <4>-<2> 

for the edges. The triangulation \i(S) then is: 



a,«i,3» = <3>-<i) 

<M(3,4)) = (4)-(3) 



H(S) = "0"®d 2 (S) = "0"® (- (1,2) + (1,3) - (2,4) + (3,4)) 
= -(0,1,2) + (0,1,3) -(0,2,4) + (0,3,4) 

Now, algebraically, (S) has the same boundary as S: 

d(n(S)) = d(- (0,1,2) + (0,1,3) -(0,2,4) + (0,3,4)) 

= -d (0,1, 2) + d (0,1, 3) -d (0,2,4) + d (0,3,4)) 
= -((1,2) -(0,2) + (0,1)) + ((1,3) -(0,3) + (0,1)) 
- ((2,4) - (0,4) + (0,2)) + ((3,4) - (0,4) + (0,3)) 
= -(1,2) + (1,3) -(2,4) + (3,4) 
= d 2 (S). 



(13) 



(14) 



(15) 



Each summand ± (0,x) in the above sum where the simplex starts with "0" has a 
simplex + (0,x) such that both cancel to zero. 



This approach, however, has a shortcoming: introducing a new vertex for each 
cell means that each sequence (c„,c„_i,. . . ,cq) of incident cells — a so-called cell tu- 
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pie — is represented by one simplex. But this can result in extremely many simplices — 
the number of cell tuples then grows more than exponentially with the dimension of 
the complex. This is also true for our approach proposed next, but at least it greatly 
reduces the number of simplices. 

5.1.1 Modifying Cohen and Hickey's Approach 

To ease the problem of the big number of simplices we use an existing bound- 
ary vertex instead of providing a new vertex for each cell. Then all boundary cells 
incident with that vertex degenerate, have a zero boundary, and can be dropped. 
We will now iteratively construct a signed simplicial decomposition as a morphism 
pi = (jl n , ... from our polytope complex ^ to a simplicial comple?<Q dimen- 
sion by dimension starting with dimension zero, This morphism is then the desired 
signed simplicial decomposition. 

We still denote the /-cells of our complex by X/. At the first step the vertices Xq 
in Co are simply labelled with integer numbers by /^o which gives a bijection 



and its linear continuation to Co — > Z[l, . . . ,#Xq\. This specifies a linear order < 
on the vertices and thereby a priori fixes an orientation of every simplex. The order 
is defined as a < b :<^=> IM)i a ) < IM)Q>) for all a,b £ Xq. We leave it open how 
permuting the labels can affect the size of the decomposition result and how we 
then can find an optimal /io- 

Each /-simplex is then an ascending sequence of i + 1 different vertices. Our 
assumptions on the polytope imply that each edge e EX\ has two different boundary 
vertices a and b and one of them is the starting vertex and the other is the ending 
vertex. This gives two matrix entries Di(e,a) = —Di(e,b) in the boundary matrix 
D\. An entry +1 indicates an ending vertex and —1 stands for the starting vertex. 
This matrix entry of the maximal vertex will become the new sign of our image edge 
going from minimal to maximal vertex. Hence pi\ is defined by 



This function inverts the edge orientation and the sign when its vertices are ordered 
against their total ordering imposed by /Xq and guarantees that the first vertex in the 
simplex (a,b) is always minimal. Note that the first two vertices are always different. 
Now this guarantee will be kept throughout every dimension. Here we denote the / th 
simplicial boundary by 5,-. 

We now show that our partially finished signed decomposition (jii , /Iq) is a mor- 
phism from the 1-skeleton (X\ ,Xq) to a simplicial complex A by showing 



1 with some salt added: As the simplices may overlap it can be disputed if this really qualifies as a 
"simplicial complex". Anyhow, it is a complex of simplices with the simplicial boundary operator. 



/io:X -> [1,...,#X ] C1N 



(16) 




(17) 
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Siojui ^jUoodi . (18) 

Let e be an edge in X\, a its starting vertex, and b be the ending vertex. This 
means d\{e) =b~a. and Di(e,b) = +1 and D\(e,a) = —I. Then for the right hand 
side of Equation ( fT8l we have: 

lM ) {di{e)) = lM ) {b-a) = Ho{b)-n {a) (19) 

In case a < bwe have for the left-hand side: 

Si(jUi(e)) = di{D { {e,b) x (po(a),^(b))) = 8i (no(a),^(b)) 

= Ho(b) - Ho(a) . (20) 

In case b < a the left hand side gives: 

Si(jUi(e)) = 5i(Di(e,a) x (^o(*),^o(a))) = 5i(- (/io(fe),/io(o)» 
= -5t({|io(ft),Mo(fl)» = -(Mo(fl)-Mo(*)) 
= ^o(A)-/io(a)- (21) 

In both cases we have /Xo(<?i (e)) = 5i (/Xi (e)) for every edge e and therefore Equa- 
tion dT8l> holds. 

The higher dimensional maps fj.2 , fl^ , . . . can now be computed iteratively each of 
them using its previously accomplished predecessor: Let /I, be given. Then we use 
it to compute jU, + 1 from /I,- and <9 I+ 1 for each i + 1 -cell c in Q+i : First take the given 
boundary <9 !+ i (c) which is a linear combination 

d f+ i(c) = ai</H h0£jfc<4 (22) 

where the a ; - are the entries Di + \{c,dj) in the boundary matrix of c?,+i. Then 
compute the image of <9,+i (c) under the morphism /x,- which gives 

mod t+ i(c) = ai/i,-(Ji)H Va k pLi{d k ). (23) 

We know that each simplex fli(dj) in above Equation d23l is a sequence (vo, ■ ■ ■ , V,) 
of vertex numbers (given by ;iio) in strictly ascending order. Therefore the first vertex 
number Vo is minimal in that simplex. We now take the minimal number a of all 
these vertices in the boundary of c and define 

a®jliO d i+ i(c) = a\ x a ® jUi(d'i) H h 0^ x a <g> ^,-(<4) (24) 

which is a linear combination of simplices (a, Vo, ■ ■ ■ , V,) which are the /-simplices 
of the triangulated boundary of c with the minimal vertex a attached at front to get 
a new i + 1 simplex. 

The operator a ® jU,-o <9, + 1 , in fact, is a morphism: 
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a ® 


// ■ O 


"1+1 V L // 


= M 


od 


+ i (c) - a ® 5,(/ii o 5,-+i(c)) 


= M 


o«9 


+i(c) - a <g> 5,-o ^ o <?i+i(c) 


= M 


od 


+ i(c)-fl®ft-i(^4+i(c)) 


= M 


od 


+ i(c)-a(g)/ii_ioO(c) 


= M 


od 


+ i(c) — a(g)0 


= M 


od 





by Eqn. © 

by /I,- being a morphism 
by d being a boundary 



(25) 



Now each boundary simplex (Vi, . . . , V„) which contains a has that vertex a at its 
front because it is the minimal vertex of all boundary simplices and Vi is the min- 
imal vertex of that simplex. Hence a = Vi and thus the simplex degenerates to 
(a, a, \>2, ■ ■ ■ , v„). But then every non-zero summand of its simplicial boundary also 
starts with (a, a, . . .). 

But such an element does not occur in 5, o /i,_i(c) because all vertices in all 
simplices in i (c) are different by precondition. By the morphism property 5; o 
Mi-l( c ) ~ Mi 0( 5;+i( c ) these elements can occur in fij o di + \ (c) only with a zero 
coefficient and therefore can be removed from fij o dj + \(c) without (algebraically) 
modifying the boundary 

Our new morphism ju,-+i does exactly this and thus guarantees the precondition 
that all vertices of a simplex are different: 

Mi+i (c) := a c ® M; (flc) ° (c) (26) 

where a c is the minimal vertex touching c and [lj is /I; after having removed all 
simplices containing a c . 

Example 5 (Triangulating our Square Again). The left hand side diagram shows a 
square S and the right hand side its triangulation. We do this step by step starting 
with Ho which enumerates the vertices in an arbitrary manner: 




Mo 



Now we can compute the edges which run from lower to higher vertex number: 
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Ml 




+<2,3> 



Remember, the the boundary of S is <?2 (S) = cd — bd — ab + ac. Then 

Hi{d 2 {S)) =H]_{cd) - H\{bd) - Hi{ab) + Hi{ac) 
= (1,3) -(2,3) -(-(2,4)) + (-(1,4)) 
= (1,3) -(2,3) + (2,4) -(1,4) (27) 

We take the minimal vertex 1, remove all simplices that start with 1, and then 
append 1 to the remaining simplices: 

H 2 (S) = - (1,2,3) + (1,2,4) (28) 



a ab 



+ s 



cd 



<2+ 2 
• 



bd 



(1,4) 



: • (1.2.4): 



(1-2) 



/ ;-(i> 2 >3); 



• 9- • 

1 (1,3) 3 



(2,3) 



Note that the 2-simplex (or triangle) (1,2,4) is oriented counter-clockwise just as 
the face S and therefore gets a positive sign whereas triangle (1,2,3) is oriented 
clockwise and accordingly gets a negative sign to indicate that is must be flipped to 
be compatible with the orientation of S. 

Finally the simplices, which are still tuples (vo, • • • , V,) of natural numbers are 
considered simplices in IR" by application of jU^ -1 , the inverse of /Xo- This gives 

(v , . . . , v,) = (jUq -1 (v ), . . . , Mo -1 (v,)). 

Now the n-d "commutative CSG" is defined as follows: If the sign of the n- 
simplex orientation matches the sign assigned to it by the morphism the (absolute) 
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simplex volume will be added to the polytope cell whereas in case that orientation 
sign and assigned sign are different the (absolute) simplex volume will be removed 
from the polytope. "Simplex orientation" of an n-simplex in IR" can easily be com- 
puted: The determinant 



det((vo,...,v„)) : 



Xi-Xq x 2 -x 

yi-yo yi-ya 

Zl—Zo Z2~Zq 



x„ — Xq 

y n -yo 

Zn —ZO 



(29) 



has n\ times the signed volume of the corresponding simplex and the sign is the 
orientation. The values x,-,y,-, . . . ,Zi are the n vertex coordinates of vertex v,. 



6 Simplex Intersection 



One application of such signed decomposition is polytope intersection or — more 
generally — polytope complex intersection: the pair-wise intersection of cells in a 
polytope complex which itself gives a new polytope complex. 

Wit those signed simplices jj. we compute the polytope complex intersection of, 
say, ^ and J(f . But how do we intersect simplices? 



6.1 Application of the Active-Set-Method 

Assume there are two simplices: an n-simplex a = (ao, . . . ,a n ) and an w-simplex 
b = (bo, . . . ,b,„). We will give a brief sketch of how the cells of the intersection 
complex of a n b can be found using a modified version of the active-set-method 
lfl6ll a well-known numerical optimisation technique. 

Each point p in a simplex (oq, . . . ,a„) has convex coordinates (ao, . . . , a„) such 

that p = OCoflo H 1- Ol n a n . In other words, these coordinates must sum up to 1 and 

neither is negative. Let q — /3o£>o H h j3„A„ be such point in the other simplex b. 

Then the pair of points (p,q) can be expressed by concatenating their barycentric 
coordinates to a vector x = (ao, . . . a„, j3o, ■ ■ ■ ,j3 m ). It represents an intersection iff 
(the square of) the distance of p and q is zero. The square of the distance can be 
computed by inserting x into 
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d 2 (x) 



( aoao 

a„ao 
~b a 



aoa„ —aobo 

a„a„ —a„bo 
-b a„ b b 



\ -b m ao ■ ■ ■ -b m a„ b m bo 



—aob m ^ 

®nbn 

bob m 

bmbm / 



(30) 



a so-called "quadratic form". The matrix entries are the dot product of the vertices 
in IR". To minimise this equation, subject to the constraints 

ao + --- + a„ = 1 = J8 H hj8 m , (31) 

the following Karush-Kuhn-Tucker-system (or KKT-system) must be solved: 



a„ao 
-boao 

-b m ao 
1 

V o 



a a n -a Q b Q ■ ■ ■ -a b m 1 \ / 0£q \ 



a n a„ -a„bo 
-boa„ b bo 

-b m a n b m bo 




1 



-a n b m 1 
bob m 1 

b,„b m 1 



A) 



A u 1 

1 Oj \A h J \lj 










(32) 



Note that the last two lines of this equation system are Equations (13 It . 

Now a vertex is called active at x if its corresponding barycentric coordinate is 
zero, and otherwise inactive. A vector is feasible iff it represents a convex combina- 
tion of the vertices, hence iff none of the barycentric coordinates is negative. A set 
of vertices is feasible if its KKT-system has a unique feasible solution. 
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To find all intersection cells of the two simplices we first compute the intersection 
vertices with the following algorithm: 

• Initialise a set S — {{a,-,£> ; } | i = Q...n,j = 0. . .ra}, the initially inactive sets, 
and initialise an empty set R (the "Result"). 

• While S is not empty do: 

- Chose one inactive set / from S and remove it from S. 

- Restrict the KKT-system of Equation f32| ) to / and solve it, say, by calling 
LAPACK's dgesv-routine J2J. "Restrict" means "remove each row and 
column x which is not in / from the equation system" (or, equivalently, set 
each such row and column x in the matrix to and then the diagonal entry 
at x to 1). 

- If the solution is unfeasible the intersection point is outside of one of the 
simplices. Then discard / and continue the loop with the next / from S. 

- Otherwise, if the Lagrangian multipliers A a and Aj, of the solution are both 
zero, the corresponding simplices intersect at the point with these barycen- 
tric coordinates. Then add / to R, because / then indicates an intersection 
vertex of the simplex intersection. 

- If the solution is feasible but at least one of the Lagrangians is not zero the 
simplices are either skew or parallel and their affine spaces do not intersect. 
Then multiply the solution with the Matrix of the unrestricted KKT-system. 
This gives (half) the gradient of the (squared) distance function at x. Then 
for each vertex V; ^ / where that gradient is negative add / + {v,-} to S. 

• When S is empty return result R. 

We can prove that this algorithm always terminates and computes every mini- 
mal inactive set for each vertex of the simplex intersection. We will not carry out 
that proof here. It consists in showing that every such inactive set has a sequence 
of vertices which can be removed one by one without rendering the remaining set 
unfeasible until it becomes one of the initially inactive sets. 

A forthcoming paper will show how we handle ill-conditioned matrices and "crit- 
ical" intersections like an intersection of a vertex with an edge of a triangle. 

Each set / in the above computed set R is a set of vertices denoting a pair of 
sides one from each simplex. This pair of sides intersects in one unique point. To 
get the higher-dimensional intersection cells we simply compute all possible unions 
of these inactive sets. 

The following diagram shows an intersection of two triangles on the left hand 
side together with the inactive sets of the intersection vertices on the right-hand 
side: 
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The vertex {1,3,4,6} indicates that it is on the intersection of edge (1,3) with edge 
(4,6), whereas {1,2,3,6} denotes the intersection of vertex 6 with triangle (1,2,3). 
Note that {1,3,4,5}U{2,3,4,5} = {1,2,3,4,5} — the union of the inactive sets of 
the two bottom vertices — gives the inactive set of the bottom line between these ver- 
tices. If we restrict the KKT-system to that set then its solution space is of dimension 
one. It contains the barycentric coordinates of all points in the affine space of that 
edge, whereas the feasible solutions denote the points on that edge itself. 



6.2 The Algebraic Simplex Intersection Boundary 

From now on we will write the sets {1,2,3,5} as "1235", i.e. a string of inte- 
gers. What follows is the set of all unions of "1345", "2345", "1346", "2356", and 
"1236": 

X = {"1345", "2345", "1346", "2356", "1236", the vertices 

"12345", "13456", "23456","12346", "12356", the edges (33) 
"123456"} the face 

This gives a topological data type (X, D) where each set is bounded by its subsets 
in X. To get a relational complex we need an algebraic boundary operator. A pos- 
sible solution could be to restrict the simplicial boundary to X. This gives indeed a 
complex boundary (see also 02)) but none which is useful for our purpose. We will 
indicate restriction to X by either <expr> \x or by striking eet sets which are not in 
X. Let us try that simple approach: Computing some edge boundary gives strange 
edges like 

d\ x ("23456") = +"345#"-"245#" + "2356" -"2346-" + "2345" 

= +"2356" + "2345". (34) 

Remember that a positive sign indicates "end point", and a negative sign is attached 
to a starting point. So we see that, instead of assigning a starting point and an end- 
point, that boundary operator assigns two endpoints "2356" and "2345" and no 
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starting point to the edge "23456". Hence, that boundary defines no orientation. 
Computing the boundary of the face and all edges gives the following result: 



1236 




On the right hand side we have inverted the edges with a negative face-edge- 
incidence to illustrate the edges cycle. This strange boundary does not define an 
orientation of the face — one part of the boundary surrounds the face in counter- 
clock-wise direction whereas the other part goes clock-wise. Both these parts start 
with a common edge that has two endpoints and end with another common edge 
with two starting points (upper left). Anyway, formally this is a valid cycle and the 
result, in fact, constitutes a complex. We suppose that we can take that boundary as 
it is for our next algorithmic step. 

As we might need an orientation for that step, however, we came up so far with 
an algorithm which iterates the dimensions starting with 1 (the edges) and attaches 
an arbitrary orientation to each cell. Because of the manifold property of each cell 
we can specify an orientation of c by merely attaching one arbitrary sign, say +1, to 
a cell-incidence R(c, d) from c to an already oriented boundary cell d. This sign then 
determines the signs of all other incidences of that cell c with its boundary cells. 

Finally the orientation of the n-cells may have to be inverted, because later we 
want that the orientation of each resulting n-cell is the product of the orientations of 
the intersecting n-simplices. By now this restricts the validity of the next algorithmic 
step to full-dimensional complexes, in other words, it might only work with n-d 
complexes in IR". 



7 Summing Simplex Intersections 

We simply sum up all these intersection complexes thereby respecting the previously 
attached signs to get the resulting intersection complex. If the intersected simplices 
have different signs then the intersection complex will be subtracted from the result- 
ing complex and otherwise added. Then the simplex boundary cells that have been 
introduced by the triangulation will cancel to zero. We will show that this at least 
works with full-dimensional complexes: 
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Let p be a point in IR" and let o = (vq, ■ ■ ■ ,v„) be an n-simplex in IR". Then p 
is almost never a boundary point and, hence, can be considered either an interior 
or an exterior point. If p is an interior point, the boundary "wraps around" p one 
or several times in some "direction". This is called the degree (or winding number) 
of a's boundary with respect to p and is nothing but the orientation of o which is 
non-zero. If p is an exterior point, however, we have winding number 0. 

But the winding number of a cell's boundary d(c) is the same as the sum of the 
signed winding numbers of the simplices' boundaries in the signed decomposition 
of c. This is also true for another triangulated n-cell z of the other complex. Let us 
denote the winding number of a cell boundary d(c) at point p by w(p,c). Then for 
the decomposition jU„(c) = 0Ci(7i H h CLk^k mto simplices a, we have 

w(p,c) = aiw(p,o l )-\ \-a k w(p,a k ). (35) 

Now p is an interior point of c n z iff both winding numbers are non-zero, or, 
equivalently, iff the product of both winding numbers is non-zero. Then we have 

p ecnz^>w(p,c)-w(p,z) ^0 (36) 

As we have a boundary for a simplex intersection o a D Oj, such that w(p,O a ) ■ 
w(p, Ob) = w(p, O a H Ob) we can compute w(p,c) ■ w(p,z) by 

w(p,c)-w(p,z) = {a,\w{p,0\) H \-a„w(p,Oi))- 

n m 

= ££a i /3 ; w(p,(7 ( n^) (37) 

i= 1.7=1 

As w(p,0) = we only have to sum up the winding numbers of non-empty sim- 
plex intersections with coefficients 0C;/3;. This is the reason for our above mentioned 
restriction on the orientation of the simplex intersection. Of course, here spatial in- 
dexing or a sweeping hyperplane approach would increase efficiency, but we will 
not particularise spatial indexing here. 

Then the refinement that stems from the signed decomposition must be undone. 
Therefore the simplex-simplex intersections will be re-composed to get the final 
result by using the inverse of the two morphisms jj.^ and jj, 

Finally, as intersecting possibly non-convex cells may give non-connected result- 
ing cells, it is also possible — but not necessary — to identify connected components 
of cell intersections if this is needed. 



8 Applications and Outlook 

Polytope complex intersection can be used to compute "topological relations", to 
overlay spatial structures like a geological formation and a drilling path or the "com- 



Polytope Complex Decomposition and Overlay 



21 



mon footprint" of two geological structures that meet at a fault. It can also be used 
to combine spatial data sets, say, a city model and a GIS data set, into one. It is 
dimension independent and, for example, can also be applied to cut temporal slices 
out of a 4D space-time complex. 

So we have shown a relational database schema for polytope complexes of any di- 
mension and an intersection algorithm for full-dimensional complexes which seems 
efficient. 

However, even if this combinatorial decomposition looks quite efficient there is 
a general complexity problem: If we combinatorially triangulate an n-dimensional 
hypercube we get n\ simplices. So the number of simplices can grow super- 
exponentially with the complexity of the triangulated polytope if there is no fixed 
dimension upper bound. In practice dimension might be bounded from above by 
5: E.g. three spatial dimensions, time, and, maybe, version history. Then our tri- 
angulation creates 5! = 120 simplices for a cube. In future we therefore want to 
study polytope complex overlay by using a signed convex cell decomposition into 
more general convex shapes. However, we suspect that the above mentioned com- 
plexity explosion with unbounded dimension is fundamental, yet another curse of 
dimensionality, and cannot be overcome by any algorithm be it smarter than the one 
presented here or not. 

Additionally, it is worthwhile to further investigate the strange results of "purely 
algebraic" boundary operators. First it should be investigated if such strange bound- 
ary could simply be tolerated when intersections are summed up. One could also 
search for an alternative to compute intersected cell orientation more deterministi- 
cally than choosing an arbitrary initial "seed" incidence sign. In particular all this 
should lead to intersect complexes of different dimension or of lower dimension 
than the embedding space. 

Also, the algorithm computes many simplex-simplex-intersections which later 
disappear at the re-composition phase and, hence, need not be computed. So it might 
save computation time if these intersections were identified in advance and not com- 
puted at all. 
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